(*  :Title:	Solution of Difference Equations via z-transforms. *)

(*  :Authors:	Wally McClure, Brian Evans, James McClellan  *)

(*
    :Summary:	This package will solve linear constant coefficient
		difference equations by means of the z-transform.
 *)

(*  :Context:	SignalProcessing`Digital`ZSolve`  *)

(*  :PackageVersion:  2.7	*)

(*
    :Copyright:	Copyright 1990-1991 by Brian L. Evans
		Georgia Tech Research Corporation

	Permission to use, copy, modify, and distribute this software
	and its documentation for any purpose and without fee is
	hereby granted, provided that the above copyright notice
	appear in all copies and that both that copyright notice and
	this permission notice appear in supporting documentation,
	and that the name of the Georgia Tech Research Corporation,
	Georgia Tech, or Georgia Institute of Technology not be used
	in advertising or publicity pertaining to distribution of the
	software without specific, written prior permission.  Georgia
	Tech makes no representations about the suitability of this
	software for any purpose.  It is provided "as is" without
	express or implied warranty.
 *)

(*
    :History:	Date started -- April 30, 1990
		Integrated into signal processing packages -- June 6, 1990
		Handle arbitrary initial conditions -- May 9, 1991

		This package started up from a special topics class.  
 *)

(*  :Keywords:	linear constant-coefficient difference equations	*)

(*
    :Source:	{Fundamentals of Digital Signal Processing}
		  by Lonnie C. Ludeman 
 *)

(*  :Warning:	*)

(*  :Mathematica Version:  1.2 or 2.0  *)

(*  :Limitation:  *)

(*  :Discussion:  *)

(*  :Functions:   PartialZ  ZSolve  *)



(*  B E G I N     P A C K A G E  *)

BeginPackage[ "SignalProcessing`Digital`ZSolve`",
	      "SignalProcessing`Digital`ZTransform`",
	      "SignalProcessing`Digital`InvZTransform`",
	      "SignalProcessing`Digital`ZSupport`",
	      "SignalProcessing`Support`TransSupport`",
	      "SignalProcessing`Support`ROC`",
	      "SignalProcessing`Support`SigProc`",
	      "SignalProcessing`Support`SupCode`" ]


If [ TrueQ[ $VersionNumber >= 2.0 ],
     Off[ General::spell ];
     Off[ General::spell1 ] ];


(*  U S A G E     I N F O R M A T I O N  *)

RightSided::usage =
	"RightSided is an option for ZSolve. \
	If True, then the solution to the difference equation will be \
	right-sided, left-sided otherwise."

ZSolve::usage =
	"ZSolve[ diffequ == drivingfun, y[n] ] solves the difference \
	equation diffequ = drivingfun, where diffequ is a linear constant \
	coefficient difference equation and drivingfun is the driving \
	function (a function of n). \
	Thus, diffequ has the form a0 y[n] + a1 y[n - 1] + .... \
	By using options, one can specify initial values; e.g., \
	ZSolve[ y[n] + 3/2 y[n-1] + 1/2 y[n-2] == (1/4)^n Step[n], \
	y[n], y[-1] -> 4, y[-2] -> 10 ]. \
	A difference equations of N terms needs N-1 initial values. \
	All unspecified initial values are considered to be zero. \
	Right-sided and left-sided solutions are possible. \
	ZSolve can justify its answers."

(*  E N D     U S A G E     I N F O R M A T I O N  *)


Begin[ "`Private`" ]


(*  Messages  *)
ZSolve::notsupported =
	"You specified initial conditions that start too far back."
ZSolve::noztransform =
	"The forcing function `` does not have a valid z-transform."
ZSolve::redundant =
	"Same initial index is assigned two different values."
ZSolve::toomany =
	"`` initial conditions given for a difference equation of order ``."
PartialZ::noztransform =
	"The difference equation `` does not have a partial z-transform."


(* introduction *)
introduction [ y_[n_], minindex_, maxindex_, diffequ_,
		minoffset_, maxoffset_, drivingfun_, initlist_ ] :=
Block [	{rstr},
	Print["Solving for ", y[n], " in the difference equation"];
	Print[ "  ", diffequ, " = ", drivingfun ];
	rstr = If [ TrueQ[rightsided],
		    StringForm["`` > ``", n, maxindex],
		    StringForm["`` < ``", n, minindex] ];
	Print[ "  such that ", rstr,
	       " and subject to the initial conditions" ];
	If [ EmptyQ[initlist],
	     Print[ "  being zero." ];
	     Print[ "  ", initlist ] ];
	Print[ " " ];
	Print[ "The first step is to substitute the initial conditions" ];
	Print[ "into the difference equation for ", n, " from ",
		minindex - maxoffset, " to ",
		maxindex - maxoffset, "." ];
	Print[ "This will give rise to ", maxindex - minindex + 1,
		" impulse function(s) which will" ];
	Print[ "be added to the right-hand side of the difference equation." ];
	Print[ "Then, the difference equation will be valid for all values" ];
	Print[ "of ", n, ", so the unilateral z-transform can used to solve" ];
	Print[ "the problem." ];
	Print[ " " ] ]


(*  homogeneous  *)
homogeneous[ drivingfun_, diffequ_, y_[n_], dialogue_, rsided_, {} ] :=
	diffequ == drivingfun

homogeneous[ drivingfun_, diffequ_, y_[n_], dialogue_, rsided_, initlist_ ] :=
	Block [	{auginitlist, homtable, indexlist, initcond, matrix, maxindex,
		 minindex, maxoffset, minoffset, newdiffequ, n0,
		 ntemp, offsetlist, order, R, step},

		(* determine the order of the difference equation *)
		offsetlist = Apply[List, diffequ] /.
			     {a_ y[n + b_] :> b, y[n + c_] :> c, y[n] -> 0};
		offsetlist = Select[offsetlist, IntegerQ];
		maxoffset = Max[offsetlist];
		minoffset = Min[offsetlist];
		order = maxoffset - minoffset;

		(* determine the maximum index of the initial conditions *)
		If [ EmptyQ[initlist],
		     minindex = maxindex = 0,

		     indexlist = ReplaceAll[ initlist,
					     (y[index_] -> value_) :> index ];
		     If [ ! SameQ[indexlist, Union[indexlist]],
			  MyMessage[ZSolve::redundant, Null] ];
		     initcond = Length[indexlist];
		     If [ initcond > order,
			  Message[ZSolve::toomany, initcond, order];
			  indexlist = Take[indexlist, -order];
			  initcond = order ];
		     maxindex = Max[indexlist];
		     minindex = Min[indexlist];
		     If [ ( maxindex - minindex ) >= order,
			  Message[ZSolve::notsupported] ] ];

		(* update dialogue with user *)
		If [ SameQ[dialogue, All],
		     introduction[ y[n], minindex, maxindex, diffequ,
			    minoffset, maxoffset, drivingfun, initlist ] ];

		(* determine effects of initial conditions *)
		newdiffequ = diffequ /. n -> ntemp;
		homtable = Table[(newdiffequ) Impulse[n - ntemp],
				 {ntemp, minindex - maxoffset,
					 maxindex - maxoffset}];
		auginitlist = Append[initlist, y[index_] :> 0];
		step = If [ TrueQ[rightsided],
			    Step[n - (maxindex + 1)],
			    Step[(mindex - 1) - n ] ];

		diffequ == drivingfun step +
				Apply[Plus, homtable /. auginitlist] ]


(*  preprocess  *)
preprocess[ left_ == right_, y_[n_], dialogue_, rs_, initlist_ ] :=
	Block [	{diffequ, drivingfun},

		diffequ = Select[left - right,
				 MatchQ[#1, a_. y[n + n0_.]]&];
		drivingfun = - Select[left - right,
				      ! MatchQ[#1, a_. y[n + n0_.]]&];

		homogeneous[drivingfun, diffequ, y[n], dialogue, rs, initlist] ]


(*  PartialZ  *)
PartialZ[ x_ + v_, rest__ ] :=
	PartialZ[ x, rest ] + PartialZ[ v, rest ]

PartialZ[ a_ x_, y_[n_], z_, YY_ ] :=
	a PartialZ[ x, y[n], z, YY ] /; FreeQ[a, n]

PartialZ[ y_[n_ + n0_.], y_[n_], z_, YY_ ] :=
	z^n0 YY /; FreeQ[n0, n]

PartialZ[ a_, y_[n_], z_, YY_ ] :=
	MyMessage[ PartialZ::noztransform, 0, a ]


(*  ZSolve  *)
ZSolve/: Options[ZSolve] := { Dialogue -> True, RightSided -> True }

ZSolve[ y_[n_ + a_.] == drivingfun_, y_[n_] ] :=
	( drivingfun /. { n -> n - a } )

ZSolve[ left_ == right_, y_[n_], init___ ] :=
	Block [	{dialogue, diffequ, drivingfun, equationlist, index,
		 initlist, options, rightsided, rminus, rplus, X, Y,
		 Ycollected, Ypartial, YY, YYstr, z, zobj, zsol},

		(* Parse options, determine the level of Dialogue requested *)
		options = ToList[init] ~Join~ Options[ZSolve];
		dialogue = Replace[Dialogue, options];
		rightsided = Replace[RightSided, options];
		initlist = Select[ options,
				   MatchQ[#1, y[index_] -> value_]& ];

		(* Collect terms, augment driving function with i.c.'s *)
		equationlist = preprocess[ left == right, y[n],
					   dialogue, rightsided, initlist ];
		diffequ = First[equationlist];
		drivingfun = Second[equationlist];

		(* Tell the user about the augmented difference equation *)
		If [ dialogue ||  SameQ[dialogue, All],
		     Print["Solving for ", y[n], " in the difference equation"];
		     Print["augmented by the initial conditions:"];
		     Print["  ", diffequ, " = ", drivingfun] ];

		(* Find complete z-transform of driving function *)
		zobj = ZTransform[ drivingfun, n, z, Dialogue -> False ];
		If [ ! SameQ[Head[zobj], ZTransData],
		     Message[ ZSolve::noztransform, drivingfun ];
		     Return[] ];
		rminus = GetRMinus[zobj];		     (* R- of X(z) *)
		rplus = GetRPlus[zobj];			     (* R+ of X(z) *)
		X = TheFunction[zobj];		     (* X(z) from z object *)

		(* Find partial z-transform of difference equation *)
		Ypartial = PartialZ[ diffequ, y[n], z, YY ];

		(* Dialogue taking z-transform of both sides to user *)
		If [ dialogue,
		     Print[ "After taking the z-transform of both sides and" ];
		     Print[ "solving for the z-transform of ", y[n], ":" ] ];
		If [ SameQ[dialogue, All],
		     Ycollected = Collect[Ypartial, YY];
		     YYstr = ToString[StringForm["Z{ `` }", y[n]]];
		     Print[ "The z-transform of the left side is:" ];
		     Print[ "  ",
		            ( Ycollected /. YY -> YYstr ) /. z -> "z" ];
		     Print[ "The z-transform of the right side is:" ];
		     Print[ "  ", X /. z -> "z" ];
		     Print[ "Solving for the unknown z-transform yields" ] ];

		(* Inverse transform complete z-transform of y, Y(z) *)
		zsol = Solve[ ( Ypartial == X ) /. z -> z^-1, YY ];
		Y = ( YY /. First[ zsol ] ) /. z -> z^-1;

		If [ dialogue || SameQ[dialogue, All],
		     Print[ "  ", Y /. z -> "z" ];
		     Print[ "Inverse transforming this gives ", y[n], ":" ] ];

		InvZTransform[{Y, rminus, rplus}, z, n, Dialogue -> False] ]


(*  E N D     P A C K A G E  *)

End[]
EndPackage[]

If [ TrueQ[ $VersionNumber >= 2.0 ],
     On[ General::spell ];
     On[ General::spell1 ] ];


(*  H E L P     I N F O R M A T I O N  *)

Combine[ SPfunctions, { ZSolve } ]
Protect[ ZSolve ]


(*  E N D I N G     M E S S A G E  *)

Print[ "ZSolve, a difference equation solver, is loaded." ]
Null
